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Abstract 

The piece- wise parabolic method (PPM) is apphed to simulations of forced isotropic 
turbulence with Mach numbers ~ 0.1 ... 1. The equation of state is dominated by the 
Fermi pressure of an electron-degenerate fluid. The dissipation in these simulations 
is of purely numerical origin. For the dimensionless mean rate of dissipation, we 
find values in agreement with known results from mostly incompressible turbulence 
simulations. The calculation of a Smagorinsky length corresponding to the rate of 
numerical dissipation supports the notion of the PPM supplying an implicit subgrid 
scale model. In the turbulence energy spectra of various flow realisations, we find 
the so-called bottleneck phenomenon, i. e., a flattening of the spectrum function 
near the wavenumber of maximal dissipation. The shape of the bottleneck peak in 
the compensated spectrum functions is comparable to what is found in turbulence 
simulations with hyperviscosity. Although the bottleneck effect reduces the range 
of nearly inertial length scales considerably, we are able to estimate the value of 
the Kolmogorov constant. For steady turbulence with a balance between energy 
injection and dissipation, it appears that C ~ 1.7. However, a smaller value is found 
in the case of transonic turbulence with a large fraction of compressive components 
in the driving force. Moreover, we discuss length scales related to the dissipation, 
in particular, an effective numerical length scale Aefj, which can be regarded as the 
characteristic smoothing length of the implicit filter associated with the PPM. 
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1 Introduction 



The paradigm of statistically stationary isotropic turbulence put forward by 
Kolmogorov is based on the balance between energy injection and dissipation 
in equilibrium, regardless of the mechanism of dissipation. This implies the 
well-known scaling relation [1] 



for the turbulence energy spectrum in the inertial subrange of wave numbers 
k. The universality of Kolmogorov scaling was confirmed in many laboratory 
experiments and three-dimensional numerical simulations [2,3,4,7,8,10,11,14]. 
The mean rate of dissipation (e) is determined by integral length and velocity 
scales, L and V, respectively, which are related to the properties of energy 
injection into the flow. Normalising (e) in terms of these characteristic scales, 
a dimensionless quantity is obtained: 



It is customary to specify the dimensionless rate of dissipation in terms of 
the similarity parameter ~ (e). The exact definition of is given in sec- 
tion 3. Numerous attempts have been made to infer the value of either 
from laboratory measurements or numerical simulations of turbulent flows. It 
appears that (e) is generally of the order unity. Particularly, a universal value 
of about 0.5 for statistically stationary isotropic turbulence at sufficiently high 
Reynolds number has emerged over the last decade [11,12,13]. However, virtu- 
ally all of the numerical estimates of have been obtained from simulations 
of incompressible turbulence so far. Only very recently, first results from sim- 
ulations of compressible flows with higher-order flnite difference methods have 
become available and, actually, are in excellent agreement with previously re- 
ported values [16]. We will extend this approach and present calculations of 
the instantaneous rate of dissipation from the energy budget in numerical sim- 
ulations of forced compressible turbulence up to Mach numbers of the order 
unity. Moreover, the whole range of developing, steady and decaying turbu- 
lence is investigated. In order to make parameter studies feasible, we chose a 
rather moderate resolution of 432^. 

The simulations we have performed also differ in other aspects from those 
cited above. Firstly, this work was motivated by the astrophysical problem of 
turbulent deflagration in degenerate stars. Thus, both the equation of state 
and the absolute values of physical quantities differ greatly from anything that 
is encountered in engineering or in the atmospheric sciences [20,21]. In terms 
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of dimensionless quantities, however, no significant deviations from turbulence 
in ideal gas are found. Apart from that, we use a high-order Godunov scheme, 
namely, the piece- wise parabolic method (PPM), which applies to fully com- 
pressible flows and, in fact, introduces numerical dissipation [5]. Of course, the 
question arises, whether the action of numerical dissipation properly mimics 
the effect of Navier-Stokes viscous dissipation in direct numerical simulations 
(DNS) or the action of subgrid scale turbulence stresses in large-eddy sim- 
ulations (LES) [9,18]. One cannot reasonably expect that the local rate of 
dissipation in the former case or the energy transfer toward unresolved scale 
in the latter case is exactly matched by the dissipative effects of the PPM. 
However, at sufficiently high resolution, a nearly inertial range of scales will 
emerge, in which the flow is dominated by the dynamics of the turbulence cas- 
cade and, hence, more or less independent of the actual dissipation mechanism. 
In addition, if statistical properties such as the value of or the Kolmogorov 
constant C were reproduced correctly in simulations of turbulence with the 
PPM, then the evolution of mean quantities and basic structural features of 
the computed flow should be in accordance with more accurate methods with 
an explicit treatment of dissipation. Indeed, there are several numerical results 
in support of this point of view [4,8,15], and we will attempt to strengthen 
the case for the PPM in this article. 

Particularly, we computed turbulence energy spectrum functions from selected 
flow realisations. The spectra for fully developed turbulence reveal a pile-up of 
kinetic energy near the wave number of maximum dissipation. This is known 
under the name bottleneck effect and was found in numerous numerical simu- 
lations [9,11,14]. The bottleneck is believed to be a genuine feature of isotropic 
turbulence and can be attributed to the partial suppression of non-linear tur- 
bulent interactions under the influence of viscous dissipation [17]. Based on the 
turbulence energy spectra, we propose the definition of a characteristic length 
scale Aeff of a presumed implicit filter which is equivalent to smoothing effect 
of the PPM. If the equations of motion are solved with spectral methods, then 
Acflf = A, where A is the length scale of the numerical grid in physical space. 
For finite volume methods, however, the smoothing of velocity fiuctuations due 
to numerical dissipation on the smallest resolved length scales / ~ A implies 
^eflf — /3A, where /3 is expected to assume a valuer larger than unity. Further- 
more, values for the dimensionless rate of dissipation 2 were calculated from 
the growth rate of the internal energy corrected for compressibility effects. 
Combining values of the numerical dissipation rate and the effective length 
scale Aeff, respectively, it was possible to calculate the Smagorinsky constant, 
assuming the the numerical dissipation would be statistically equivalent to the 
action of a Smagorinsky subgrid scale model. 

In section 2, we will outline the forcing scheme, the simulation parameters 
and the equation of state in our simulations. The evaluation of the rate of 
dissipation is presented in section 3, the computation of turbulence energy 
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spectra in section 4 and, finally, dissipation length scales are discussed in 
section 5. 



2 Forced isotropic turbulence 



The dynamics of forced turbulence in compressible fluids is determined by the 
following set of conservation laws in combination with the equation of state 
[22,23]: 

9 d 

Ft^ + = P^''^'- 

The field f{x,t) is called the driving force and may be any mechanical force 
acting upon the fluid and thereby supplying energy to the flow. The term aik,k 
in the momentum equation is the viscous force per unit volume. In the case of 
Navier-Stokes turbulence, the viscous stress tensor is proportional to the 
local strain of the velocity field: 

aik = "^P^S-k = 2pu (^Sik - ^Vjj6ik^ , (6) 

where Sik — ^{vi^k + Vk,i)- The coefficient u is the microscopic viscosity of the 
fluid. 



According to the Landau criterion [22], the range of dynamically relevant 
length scales in a turbulent flow is of the order L/rjK ^ Re^''^, where L is the 
integral length scale associated with the spatial variation of the driving force, 
and rjK, which is called the Kolmogorov scale, is a characteristic length scale 
associated with viscous dissipation. If the Reynolds number Re = LV/ v for a 
flow of characteristic velocity V is large enough, the range of dynamically rele- 
vant length scales will become computationally intractable. Then the common 
approach is to run a LES which resolves only the largest scales and invokes 
some model for turbulence on smaller scales. On the other hand, it has been 
suggested to let merely the dissipative effect of a finite- volume scheme smooth 
out the flow on length scales just above the numerical cutoff and dispose of ac- 
cumulating kinetic energy [18] . In particular, this assumption was put forward 
and exhaustively tested for the piece-wise parabolic method (PPM) [5,7,8,9]. 

We chose to follow this approach and adopted the PPM for simulations of 
forced isotropic turbulence. In order to produce turbulence, a random driv- 



4 



ing force is applied [24]. The evolution of the force field is determined by a 
three-component Ornstein- Uhlenbeck process in spectral space [25], which cor- 
responds to the following stochastic differential equation (SDE) for the Fourier 

transform f{k, t): 

5{k-kjimmik)-dm. (7) 

/ 

The discretisation in physical space induces a discrete spectrum of modes asso- 
ciated with the wave vectors kji„i. Gaussian random increments are introduced 
by the three-component Wiener process and are projected by means of the 
symmetric operator 

mdk) = C^.^(^) + (1 - CMik) = CSi, + (1 - 20^. (8) 

For the spectral weight C = 1; the physical force field becomes purely solenoidal, 
i. e., divergence-free. Choosing C < 1; dilatational components are generated 
which compress or rarify the fluid. The spectral profile of the driving force 
is determined by the variance a'^(k). We use a symmetric quadratic function 
centred at the wave number ko- For k > 2kQ, all modes of the force vanish 
identically. The integral length scale L is defined by L = 27r/A;o. Once a{k) 
is normalised, the asymptotic RMS amplitude of the stochastic driving force 
depends on Fq and the weight ( only: 

f,^ ~ (1 - 2C + 3C')Fo. (9) 

The linear drift term in the SDE (7) causes any information about the initial 
conditions to decay exponentially over the characteristic time T [25]. Con- 
sequently, the computed fiow will evolve towards a statistically stationary 
state which becomes asymptotically independent of the initial state. In this 
regime, the properties of the flow are determined by the three parameters L, 
Fq and C- The force magnitude Fq is conveniently expressed as Fq = V^/L, 
where the characteristic velocity V is about the root mean square velocity 
in the fully developed flow, and for the integral time scale we naturally set 
T=(L/Fo)V2 = L/y. 

Civen the astrophysical background of our work, we employed the equation of 
state (EOS) of a degenerate electron gas. Degeneracy is a quantum mechan- 
ical phenomenon in the limit of vanishing temperature which, for example, 
is encountered in so-called white dwarfs. This kind of stellar remnant em- 
anates from the burnout of stars comparable in mass to our Sun [20] . In white 
dwarf matter, the electrons are not bound to nuclei. While the latter are non- 
degenerate and obey the ideal gas law, each electron assumes more or less 
its lowest possible energy state. Thus, the electrons follow the Fermi-Dirac 
statistics. One can show that the resulting degeneracy pressure of the electron 
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gas dominates over the residual thermal pressure of all particles in the inte- 
rior of a white dwarf [27]. In this case, the equation of state is approximately 
given by P oc p^^^ independent of the temperature, provided that the mass 
density is less than the critical density Pc ~ rrip/Xl pa 2.9 • lO^^kgm"^, where 
nip is the mass of the proton and Ae = h/rrieC is the Compton wave length 
of the electron. At higher densities, the electrons become relativistic and the 
isentropic exponent approaches 4/3 [27]. 

In fact, a realistic EOS which accounts for contributions from degenerate elec- 
trons, non-degenerate nuclei, photons and pair production over a wide range 
of densities and temperatures was used in the simulations we performed. This 
EOS cannot be expressed in analytical form and must be evaluated numeri- 
cally [28]. The following initial conditions were chosen in all but one case: 

v{x, 0) = 0, (10) 
p(x, 0) = po = 0.02-| 5.805 • 10^ kg m'^, (11) 

2 

T{x, 0) = To = 0.001-^ « 5.929 • 10^ K. (12) 

Kb 

For the simulation with the lowest Mach number, the initial mass density is 
larger by a factor 2.5-10^ than the above value. In addition to po, the spectral 
weight C of the driving force, the characteristic velocity V and the integral 
times scale T as well as supplementary simulation parameters are listed in 
table 1. The spectrum of the driving force is determined by the wave number 
ko = 271 /L = 27ia/X, where a = 3 is the ratio of the domain size X to the 
integral length L. We used a grid of 432^ cells for each simulation, and an 
integral length L — 1.44 • 10^ m. This corresponds to a numerical resolution of 
A = 10.0 m. 
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Table 1 

Chosen values for basic simulation parameters, the time of the onset of decay t^/T 
and the end of each simulation t{/T, as well as the total number of time steps N/\t- 

From the hydro dynamical point of view, the EOS for degenerate matter is 
remarkable, because the pressure is virtually independent of the temperature. 
In consequence, thermodynamics and hydrodynamics are decoupled. In the 
case of an ideal gas, on the other hand, both the temperature and the pressure 
are gradually rising in a turbulent flow due to the heat generated by kinetic 
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Fig. 1. Evolution of dimensionless mean quantities for a simulation of forced isotropic 
turbulence in three dimensions. The characteristic Mach number is V/cq = 0.66 
and the spectral weight of the driving force C = 0.75. The panels show the RMS 
momentum and specific force, the mean total and internal energy, the RMS vorticity 
and divergence as well as the averaged rates of energy production and dissipation 
as functions of the normalised time t = t/T. 

energy dissipation. Strictly speaking, turbulence can never be stationary in an 
ideal gas, because of the changing relation between fluctuations of the density, 
pressure and velocity, respectively, with varying temperature. Of course, there 
are limitations in the case of degenerate matter as well, because the rising 
internal energy will eventually break the degeneracy of the electron gas. This 
happens once the temperature becomes comparable to the Fermi temperature, 
which is > 10^ K. In the case of the simulations discussed subsequently, the 
total energy injected into the system is a significant fraction of the initial 
internal energy, but degeneracy is maintained over several integral time scales. 



The global statistics of the simulation with V ~ 0.66co, where Cq is the initial 
speed of sound, and C, = 3/4 is shown in figure 1. In this case, the major 
component of the driving force is solenoidal, but there is a significant dilata- 
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Fig. 2. Statistics for a simulation with V jc^ = 0.66 and Q = 0.75. 

tional contribution as well. The plot of the RMS velocity shown in the top left 
panel illustrates the different regimes of the flow evolution. In the beginning, 
fluid is gradually set into motion. In the course of production, the intensity of 
turbulence grows exponentially, as one can see from the evolution of the RMS 
vorticity = |V x ul in the bottom panel on the left. At time t = t/T ^ 2, 
the increase of the kinetic energy stagnates. From t ^ 3 onwards, the mean 
kinetic energy remains roughly constant, and there is an approximate balance 
between energy injection and dissipation (panels on the right of figure 1). 
Moreover, the graph of the mean vorticity flattens for t > 2. Altogether, we 
take the evolution of the statistical quantities as phenomenological indication 
of fully developed turbulence at time t > 2, which asymptotically approaches 
a statistically homogeneous and stationary state. The conspicuous filaments 
of intense vorticity, which are the hallmark of developed turbulence, can be 
seen in a three-dimensional visualisation of the flow at time t = 4.0 in the 
figure 7. The intermittency of turbulence becomes manifest in the spacious 
voids between the vortex filaments. At time t = t^ = 5.0, the dropping of the 
random diffusion term in the SDE (7) initiates the exponential decay of the 
force field, and the energy contents of the flow is subsequently dissipated. 
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Quite a different behaviour emerges in the case V/cq = 1.39 and ( = 1/5. 
As one can see from the plots in figure 2, the rise of the RMS velocity is 
less steep. Actually, the evolution of the RMS structural invariants (l^*!^)^/^, 
(a;^)^/^ and {cP)^/'^, which is shown in the left bottom panels, suggest that 
there is an initial phase which is dominated by shocks rather than eddies. 
This is reflected in pronounced oscillations in the rate of change of kinetic and 
internal energy, respectively, due to compression and rarefaction effects. At 
i ~ 2, vorticity begins to dominate over the divergence. Apparently, there is 
a transition from the shock-dominated phase toward a regime, in which the 
flow becomes increasingly solenoidal, despite of the mostly dilatational force 
acting upon the fluid. However, since the decay regime was initiated at the 
time t = 5.0 too, a statistically stationary flow was not established in this 
simulation. 



Furthermore, we computed a couple of mean structural invariants which cor- 
respond to third-order statistical moments of the rate of strain tensor Sij. 
Normalisation with respect to the mean rate of strain l^"! = {2SijSijy^^ gives 
rise to the definition of the following parameters: 



{p\sr\s\) 

Po(|5P)3/2' 

{pS,jSji,Ski) 



(13) 

skew(5fe). (14) 



The first parameter, a, is closely related to the mean rate of dissipation in the 
Smagorinsky model. We will comment on this point further in the following 
section. The parameter b, on the other hand, is proportional to the skewness 
of the rate of strain tensor. Numerical values of a and b evaluated from flow 
realisations in three simulations at certain instants of time are summarised 
in table 2. Thereby, a remarkable similarity is revealed. In particular, the 
rate-of-strain skewness of about —0.3 agrees with known results for isotropic 
turbulence [26] . We take this as a further indication that the flow reahsations 
in our simulations are in good approximation both statistically stationary and 
isotropic after a few integral time scales have elapsed. 
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Table 2 

Normalised third-order moments of the rate of strain defined in equations (13) 
and (14) and the equivalent Smagorinsky length in units of A and Aeff , respectively. 
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3 The rate of dissipation 



The details of the numerical dissipation produced by the PPM are basically 
unknown, but it is possible to infer the mean rate of dissipation from the 
globally averaged energy conservation laws [23]. In fact, there are two distinct 
mechanisms of dissipation in a compressible fluid. On the one hand, pressure- 
dilatation accounts for the conversion of internal energy into mechanical en- 
ergy and vice versa, as fluid, respectively, expands or contracts. Although 
pressure-dilatation might locally produce mechanical work, it is effectively a 
dissipative, irreversible process, and the net rate of heat production is given 
by —{Pd), where d = Su is the divergence of the velocity field. If pronounced 
shocks are present, however, there might be transient phases in which {Pd) 
becomes positive. This can be seen, for example, in the right bottom panel of 
figure 2, where the time derivative of the internal energy exhibits local minima 
corresponding to the production of mechanical energy by pressure-dilatation. 
On the other hand, the change of internal energy caused by viscous dissipation 
is strictly negative. 

For periodic boundary conditions, the total flux through the boundary surfaces 
cancels out, and averaging of the internal energy equation thus yields 

d 

Poenura = -^{Eint) + {Pd). (15) 

Since we do not employ an explicit viscosity term in the equation of motion, 
enum is indeed the mean rate of dissipation due to numerical effects. The corre- 
sponding dimcnsionlcss dissipation rate inum is defined by equation (2). Several 
representative values of enum and the ratio — (P(i)/enum at certain instants of 
time are listed in the tables 3, 4 and 5. 

In the literature, it is common to normalise the rate of dissipation in terms 
of the RMS velocity fiuctuation v' and an integral length scale L, which is 
defined by the transversal turbulence energy spectrum E{k)^: 

2v'^ Jo k ~ 4v'^ ^ kn 

The average kinetic energy $^ contained in Fourier modes of wave number 
kn will be defined in the next section. For isotropic turbulence, v'"^ — |(ekin), 
where Ckin = |f ^ is the specific kinetic energy. The parameter specifying 
the dimensionless rate of dissipation is then given by 

Ce=^(enum). (17) 

Using the instantaneous values of the rate of dissipation and the RMS velocity 
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fluctuations, we estimated for different flow realisations. The accuracy is 
limited by the calculation of L, which is mostly determined by the smallest 
wave numbers, and only few discrete cells are available for A; ~ 1 in Fourier 
space. The results are listed in the tables 3, 4 and 5. Generally, Cg is small 
compared told unity in the production phase. For steady turbulence, values in 
the range 0.4 . . . 0.5 are found, which are close to the time-averaged asymptote 
~ 0.5 in simulations of incompressible turbulence [11,12,13]. 

In the introduction, we mentioned that simulations of Eulerian fluids with the 
PPM, to a certain degree, are equivalent to a LES with explicit modelling of 
the energy transfer toward unresolved scales in place of numerical dissipation. 

The consistency of this assumption in terms of statistical quantities can be 
corroborated for a simple algebraic subgrid scale model such as the Smagorin- 
sky model. The closure for the rate of dissipation in the Smagorinsky model 
is given by 

(0 = a4(|^r)^/^ (18) 

where the cocfiicicnt a is defined by equation 13, \S\ is the rate of strain and 
£s is a characteristic length of the model. Invoking the condition that (enum) 
equals the rate of dissipation predicted by the Smagorinsky model for the 
given average rate of strain in a particular flow realisation, the Smagorinsky 
length £s can be calculated. The obtained values in units of A for developed 
turbulence of three different characteristic Mach numbers are listed in table 2. 
Moreover, we normalised £s in terms of the effective length scale Agg. The 
procedure for the calculation of Aeff is introduced in section 5. As we shall 
argue, Agff is the appropriate cutoff scale for the PPM. Indeed, the results 
listed in table 2 verify in each case that £s/^cfr is about 0.16, which is just 
the Smagorinsky constant Cq obtained from analytical considerations [25]. 
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Table 3 

Mean energy, rate of dissipation and characteristic length scales for a DNS with 
V/cQ = 0.42 and C = 1-0. 
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Table 4 

Mean energy, rate of dissipation and characteristic length scales for a DNS with 
V/cq = 0.66 and C = 0.75. 
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Table 5 

Mean energy, rate of dissipation and characteristic length scales for a DNS with 
V/cq = 1.39 and C = 0.2. 

4 Turbulence energy spectra 



In the case of isotropic turbulence, the sum over all squared Fourier modes of 
the velocity field can be expressed as an integral over a function of the wave 
number only: 

77^ E i;{^3im{t) ■ v.Ut)) = / d{ak)E{ak, t) (19) 

jlm ^ 

Here the wave number is written in dimensionless form as A; = Lk/2iT. The 
function E{k,t) = {aL/27i)V^E{ak,t) is called the energy spectrum function. 
For a discrete spectrum of modes, E{ak,t) is a generalised function which is 
defined by: 

Eiak, = 7^ E \i^Mt) ■ ^Ut))SH - if + + mY']- (20) 

jlm 

For the numerical evaluation, we define an approximation to E{ak, t) by sum- 
ming over wave number bins [kn~i/2, kn+i/2]- This yields the discrete set of 
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numbers 

EM = (21) 

where kn = \{kn-i/2 + ^n+1/2), and $n(^) is the average kinetic energy of all 
modes in the corresponding wave number bin: 

^n{t) = Y E E \%im{t) ■ ^vUt)4?P^m^- (22) 

n— 1/2 — — n+1/2 

Here is the statistical weight of the nth wave number bin, i. e., the number 
of cells in Fourier space with wavenumber kjim £ [^n-1/2, ^n+i/2]- 

A number of spectra were computed from flow realisations at selected times. 
Some of these spectra are shown in Figures 3, 4 and 5. For the calculation of 
the discrete spectrum function as given by definition (21), a coarse equidistant 
wave number mesh was used in the energy containing range, < A; < 2, and 
a logarithmic mesh with narrow bins for the larger wave numbers. The panels 
on the top of figure 3 show, from left to right, the discrete energy spectrum 
function En{t) at representative stages for the simulation with C, — 1.0 and 
characteristic Mach number V/cq ~ 0.42. The panel on the very left shows 
a spectrum in the production phase. The middle panel corresponds to de- 
veloped turbulence, and the panel on the very right shows a spectrum for 
decaying turbulence. One can clearly see the larger energy contents at high 
wave numbers in the case of developed turbulence and the lower overall energy 
budget in the decay regime. Also plotted are the longitudinal and transversal 
spectrum functions, E\(f) and E^{t), respectively. El[{f) corresponds to the 
dilatational (rotation-free) components of the fiow and E^(t) to the solenoidal 
(divergence-free) components. From the plots in figure 3, it becomes apparent 
that the longitudinal energy fraction at higher wave numbers decays faster in 
comparison to the transversal fraction. The energy in the longitudinal com- 
ponents at small wave numbers, on the other hand, remains almost constant. 

From the Kolmogorov spectrum function (1), it follows that {e)~'^/^k^/^E-^[k) 
is approximately constant in the inertial subrange of wave numbers. This sug- 
gests the definition of a compensated spectrum function, 

—2/3 

{aKfl^E^{t), (23) 

as an indicator of Kolmogorov scaling. Note that only the transversal part 
of the energy spectrum is compensated, because it is the incompressible frac- 
tion of turbulence energy which a priori fulfils Kolmogorov scaling. For the 
calculation of ^;|;(t), the values of (cnum) listed in the tables 3, 4 and 5 were 
substituted for the dimensionless mean rate of dissipation in equation (23). 
The resulting plots of the compensated spectrum functions for C, = 1.0 are 
shown in the bottom panels of figure 3. In particular, the graph of ^ni^) 



a 
2^ 



m) 
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t/T =1.50 t/T = 3.00 t/T = 6.00 
1 t ' ' ' ' 'i ^ ^ ' ' 'I ^ 




(L/27T)k (L/2Ti)k (L/2Tr)k 

Fig. 3. Turbulence energy spectrum functions for V/cq = 0.42 and = 1.0. In the 
bottom panels, the compensated spectrum functions corresponding to the transver- 
sal component are plotted. 

t = 3.0 verifies the existence of a narrow window of wave numbers in the 
vicinity of k = 3.0, in which Kolmogorov scahng with C ~ 1.7 appHes within 
the bounds set by the numerical uncertainty. It is noteworthy that this value 
is very close to results obtained from other numerical simulations, especially, 
those with higher resolution [4,10,11]. However, the Kolmogorov constant in- 
ferred from our simulations might be systematically too large, because of the 
lack of isotropy for wave numbers close to the energy-containing subrange [4]. 
Moreover, the portion of the compensated spectrum function \l/;J;(t = 3T) in 
the range of dimensionless wave numbers between 2 and about 4 appears to be 
consistent with a modification of the Kolmogorov exponent by —0.1 , which 
has recently been inferred from a DNS with extremely high resolution [11]. 

Basically the same results hold for the simulation with higher characteristic 
Mach number, V/cq ~ 0.66, and ( = 0.75. The sample spectra shown in 
figure 4 are very similar in shape compared to those in figure 3, except for the 
smaller gap between the total and the longitudinal energy spectrum functions. 
As for the case ( = 0.2 and V/cq ~ 1.39, differences can evidently be seen in 
figure 5. At early time, the flow is dominated by large shock waves and most of 
the energy is contained in longitudinal modes. This is illustrated by the plot of 
the energy spectrum functions at time t = 2.0. The scaling deviates markedly 
from the Kolmogorov law. In particular, the longitudinal component 
clearly obeys a power law with exponent —2 and the transversal component is 
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(L/27i)k 



(L/27T)k 



Fig. 4. Turbulence energy spectra for F/cq = 0.66 and ( = 0.75. 

t/T = 2.00 t/T = 5.00 t/T = 10.00 




1 10 
(L/27i)k 



1 10 

(L/2TT)k 



Fig. 5. Turbulence energy spectra for V/cq = 1.39 and C = 0.2. 

falling off even steeper. From statistics of the velocity field shown in figure 2, 
a transition to the regime of solenoidal turbulent flow around t ~ 3.0 can be 
discerned [6] . In this regime, the small-scale dynamics is eventually dominated 
by turbulent vortices as in the case of subsonic turbulence, and the transversal 
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component of the energy spectrum function is more or less Kolmogorovian. 
However, we have a somewhat smaller value of the Kolmogorov constant, 
C ~ 1.3, at time i = 5.0 (see the middle panels in figure 5). 

For developed turbulence, there is a pronounced maximum of the compensated 
spectrum function at A; ^ 15 corresponding to a flattening of the energy spec- 
trum in comparison to the Kolmogorov law. This so-called bottleneck effect 
has actually been observed in various numerical simulations [3,9,14,15]. The 
anomalous bottleneck scaling is attributed to dynamical peculiarities on scales 
which are significantly influenced by dissipation. For example, experimental 
data indicate a power law behaviour of the energy spectrum function in 
the vicinity of the wave number of maximum dissipation [2]. A theoretical 
explanation of the bottleneck effect on grounds of the non-linear turbulent 
transfer was suggested too [17]. The peak of ^n(^) close to the wave number 
/(^ 15 is more pronounced than in simulations of incompressible turbulence 
with spectral methods. For instance, the peak value in the middle panels of the 
figures 3 and 4 is about 2.9, whereas the peaks of the compensated spectra in 
[3,11] are not much higher than 2.0. The width of the bottleneck, on the other 
hand, is about one decade in wave number space in all cases. Remarkably, 
the properties of the bottleneck we observe are very similar to what is found 
in turbulence simulations with hyperviscosity [15]. Thus, it appears that the 
numerical viscosity of the PPM acts like a hyperviscosity, at least as far as 
the spectral distribution of turbulence energy is concerned. 



5 Dissipation length scales 

One can regard the action of numerical dissipation as being equivalent to an 
implicit filter smoothing the flow on a certain length scale, say, Aefr. Fourier 
modes of wave number larger than tt/ A^s are suppressed. It has already been 
pointed out in section 3 that the ratio j3 = Aeg/A is of particular interest 
for subgrid scale models, which depend on the reliable specification of the 
cutoff length. In combination with a finite-volume method, the cutoff length 
may very well be different form the grid resolution A. Although the common 
point of view holds that it is not sensible to apply a subgrid scale model in 
combination with dissipative schemes such as the PPM, subgrid scale models 
are used in some applications without active coupling to the resolved flow. 
A particular example is the modelling of turbulent combustion, where the 
subgrid scale energy is treated as a passive scalar and is used as input for a 
turbulent flame speed model [19,20,21]. 

The determination of the effective length scale Acs for the PPM makes use of 
the notion of a characteristic fllter scale [29]. For a one-dimensional fllter of 
explicitly known functional form, the characteristic length can be calculated 
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from the second moment of the Fourier transform of the filter kernel, the so- 
called transfer function. However, since the implicit filter associated with the 
PPM is not explicitly known, we have to resort to the numerically computed 
energy spectra. Using Kolmogorov's law as reference spectrum function, we 
define an effective transfer function 

where E^{k,t) = Ce^^^{t)k~^/^ . The second moment of the squared transfer 
function is given by^ 

M^'\Gl^]= r k'Gl^{Kt)dk. (25) 
Jo 

Of course, the numerically computed spectrum function E{k, t) does not con- 
form with the Kolmogorov law at wave numbers comparable to ko = ^n/L, 
because of energy injection on the largest length scales. However, the contri- 
bution of these wave numbers to the above integral is small due to the factor 
k^. Thus, we will ignore this error. 

Discretising the transfer function in the fashion outlined in section 4 and 
cutting off at the wavelength tt/A, the following approximation to the second 
moment is obtained: 

The second moment of the filter transfer function has the dimension of inverse 
length cubed. Hence, a length scale is given by (M*^^^ [G^])"^/^ and is custom- 
arily normalised with respect to the second moment of the sharp cutoff filter. 
The transfer function of this filter is given by G/\{k) = 6{7!-/A — A;), and the 
second moment is M^'^^[Gl] = {ir/Af/S. Setting M(2)[G'2^] = (7r/Aefr)V3, 
the effective filter scaling factor of the PPM is therefore estimated to be 



_ Aeff _ iV 

^ A 2 



-1/3 

(27) 



where Uc = max{n|A;„ < kc}- 



We calculated j3 for each simulation at several instants of time. The obtained 
values listed in the tables 3, 4 and 5 demonstrate that Aeff ~ 1.6 A for fully 
developed turbulence. Only for the lowest Mach number, V/co pa 0.084, the 



^ Usually, filter length scales are defined by the second moment of G rather than 
G^. Computationally, however, it is preferable to use the square of the transfer 
function. 
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scaling factor /3 1.8 was found in the stationary regime. Consequently, 
the effect of numerical dissipation appears to become more pronounced for 
decreasing Mach number. Of course, the dependence of j3 on the numerical 
resolution should be investigated too. Changes can be expected towards lower 
resolution, because the energy-containing and the dissipation subrange will 
increasingly overlap. At higher resolution, on the other hand, (3 should asymp- 
totically approach a value independent of A^. Unfortunately, the validation of 
this conjecture would require an undue amount of computational resources. 

In fact, Aeff is much smaller than the length scale of maximum dissipation, Zp, 
which is given by the maximum of klEn{t) oc k]/-^^!n{t). For fully developed 
turbulence, the peak of dissipation was found to be located close to the second 
maximum of ^'n(^), with a typical value /p ~ 0.065L pa 9 A (see tables 3, 4 
and 5). Thus, Zp can be considered as the characteristic length scale of the 
bottleneck effect. The morphology of the flow on length scales Z ^ Zp is illus- 
trated in figure 6. Shown are the isosurfaces of vorticity which correspond to 
the 97 % level of the probability distribution function. In order to suppress ve- 
locity fluctuations on scales smaller than /p, a Gaussian filter of characteristic 
length lOAeff ^21^ was applied to the flow reahsation at time i — 4.0. The 
characteristic Mach number is V/cq — 0.66 and the spectral weight C = 3/4. 
The smoothing of the velocity field over a length 2Aeff ~ /p/3 results in the 
visualisation shown in figure 7. Although the length scales smaller than Zp are 
subject to significant numerical dissipation, there is nevertheless a great wealth 
of substructure on these scales. It becomes clear that the tube-like structures 
in figure 6 are actually concentrations of smaller vortices. The magnitude of 
the flow velocity is colour-coded in the plots. It appears that a significant 
fraction of high velocity fiuctuations is present on the length scales / < /p in 
the bottleneck subrange of the turbulence energy spectrum. 



6 Conclusion 

The computation of turbulence energy spectrum functions reveals several im- 
portant properties of turbulence in numerical simulations with the PPM: 
Firstly, the range of length scales approximately satisfying Kolmogorov scal- 
ing, secondly, characteristic scales associated with numerical dissipation, and 
thirdly, the so-called bottleneck effect, i. e., an excess of kinetic energy in 
modes affected by dissipation. 

For the simulations of forced isotropic turbulence discussed in this paper, there 

is only a marginal inertial subrange. Nevertheless, we were able to estimate 
the Kolmogorov constant for samples of the fiow in the nearly statistically 
stationary regime to be C ~ 1.7, which agrees very well with published results 
from more elaborate simulations. Only in the case of a simulation with Mach 
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This figure is not available in the e-print version. 

Fig. 6. Isosurfaces of vorticity for the same flow realisation as in figure 6, however, 
with the velocity field smoothed on the length scale lOAgg = 16. 2A, which is about 
the length scale of peak dissipation. 

number of order unity and forcing dominated by dilatational components, 
a significantly smaller value of C was found. On the basis of the notion of 
an implicit filter, we define the length scale A^s specifying the characteristic 
length of smoothing due to the numerical scheme. The ratio (3 = Aeff/A 
appears to be fairly universal for statistically stationary turbulence. In the case 
of the PPM, /9 fa 1.6. Furthermore, it was demonstrated that the computed 
spectra exhibit a pronounced bottleneck effect. The corresponding peak is 
higher than in spectra obtained from simulations with spectral methods. The 
magnitude of the bottleneck effect appears to be similar to what is obtained 
in simulations with hyperviscosity. 

The parameter specifying the dimensionless rate of numerical dissipation 
assumes values of about 0.5 for developed turbulence. This is consistent with 
the time-averaged asymptotic value in the limit of large Reynolds numbers 
calculated from simulations of incompressible turbulence. Prom the mean rate 
of numerical dissipation, one can also compute an equivalent Smagorinsky 
length. The ratio of this length to Acs is very close to the analytical value of 
the Smagorinsky constant Cs. This result supports the presumed statistical 
equivalence of numerical dissipation in simulations with the PPM and the 
subgrid scale energy transfer in proper LES. 
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This figure is not available in the e-print version. 

Fig. 7. Isosurfaccs of vorticity for developed turbulence in the simulation with 
F/co = 0.66 and ^ = 0.75. The velocity field was smoothed with a Gaussian fil- 
ter of characteristic length 2Aeff = 3.24A. The colour coding corresponds to the 
dimensionless velocity vjV . 

In essence, the presented results suggest that the apphcation of the PPM in 
fluid dynamical simulations yields a fair numerical representation of turbulent 
flows, for which crucial statistical parameters are in good agreement with the 
results obtained with more accurate methods. 
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